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ABSTRACT: The lowest-lying glueballs are investigated from Nş = 2 QCD study on 
aniostropic lattices. Only the gluonic operators built from Wilson loops are involved in 
calculating the corresponding correlation functions. In the tensor channel, we obtain the 
ground state mass to be 2.367(35) GeV and 2.380(61) GeV at Mmr ~ 938 MeV and 650 MeV, 
respectively. In the pseudoscalar channel, when we use the gluonic operator whose contin- 
uum limit has the form of €;;,7'rB;Dj;B,, we obtain the ground state mass to be 2.559(50) 
GeV and 2.605(52) GeV at the two pion masses. These results are in agreement with the 
corresponding glueball masses in the quenched approximation and show little dependence 
of pion masses. In contrast, if we use the topological charge density as field operators 
for the pseudoscalar, the masses of the lowest state are much lighter (around 1GeV) and 
compatible with the expected masses of the flavor singlet gg meson. This provides a strong 
hint that the operator ¢;;,7'rB;D; B, and the topological charge density couple very differ- 
ently to the glueball states and gg mesons. In the scalar channel, the ground state masses 
extracted from the correlation functions of gluonic operators are determined to be around 
1.4-1.5 GeV, which are close to the ground state masses from the correlation functions of 
the quark bilinear operators. In all cases, the mixing between glueballs and conventional 
mesons remains to be investigated in the future. 
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1 Introduction 


Due to the self-interactions among gluons, Quantum Chromodynamics (QCD) admits the 
existence of a new type of hadrons made up of gluons, usually called glueballs. Glueballs 
are of great physical interests since they are distinct from the conventional gq mesons 
described in the constituent quark model. Glueballs have been intensively studied by QCD 
inspired phenomenological models [1-6] and also by lattice QCD. Early lattice QCD studies 
in the quenched approximation show that the lowest pure gauge glueballs are the scalar, the 
tensor, and the pseudoscalar glueballs, with masses of 1.5-1.7 GeV, 2.2-2.4 GeV, and 2.6 
GeV, respectively [7-9]. For the masses of the scalar and tensor glueballs, some preliminary 
unquenched lattice studies have given compatible results [10-13]. However, for the mass 
of the pseudoscalar glueball, a consensus has not been reached. Experimentally, there are 
several candidates for the scalar glueball, such as fo(1370), fo(1500), fo(1710), however, 
none of them has been unambiguously identified as a glueball state. On the other hand, 
J/w radiative decays are usually regarded as an ideal hunting ground for glueballs. A few 
lattice studies have been devoted to the calculation of the radiative production rate of 
the scalar and the tensor glueballs in the quenched approximation [14, 15]. The predicted 
production rate of the scalar glueball is consistent with that of fo(1710), and supports 
fo(1710) to be either a good candidate for the scalar glueball or dominated by a glueball 
component. The predicted production rate of the tensor glueball is roughly 1%. It is 
interesting to note that the BESIII Collaboration find that the tensor meson f2(2340) has 
large branching fractions in the processes J/y —> ynn [16] and J/y > yog [17]. 

Even though the quenched lattice QCD studies have provide some information on the 
existence of glueballs, it is highly desired that full lattice QCD studies can be performed 
in the glueball sector. In this work, we attempt to address this question using the Ns = 2 


Table 1. Parameters of configurations 


6 DxT £ as Mr Neonf 
2.5 128x128 5 0.114fm 650 MeV 4800 
2.5 123x128 5 0.118fm 938 MeV 10400 


gauge field configurations that we generated on anisotropic lattices. The advantage of the 
use of anisotropic lattice is two-folds: on the one hand, a large statistics can be obtained 
by a relatively low cost of computational resources, on the other hand, the finer lattice 
spacing in the temporal direction can provide a better resolution for the signals of physical 
states. As the first step, we will focus on the lowest-lying glueball states, such as the scalar, 
the tensor and the pseudoscalar states. We will also explore the quark mass dependence 
of the masses of these states by using gauge configurations with different bare quark mass 
parameters. Special attention is paid in the pseudoscalar channel, since the results of the 
previous full-QCD studies in this channel are still ambiguous. On the other hand, owing 
to the U4 (1) anomaly, gluons can couple strongly to the flavor singlet pseudoscalar meson 
(n in the Np = 2+ 1 case). A recent Nf = 2+ 1 lattice study showed that y’ could 
be extracted from the negative tail of the correlation function of the topological charge 
density operator [18]. In contrast, a similar study in the quenched approximation found a 
pseudoscalar with a mass compatible with that in the pure gauge theory [19]. Therefore, 
we would like to check if a pseudoscalar glueball state also exists apart from the 7’ state 
in the presence of sea quarks. 

This paper is organized as follows: Section 2 contains a brief description for the gen- 
eration of gauge field configurations. Section 3 presents the calculation details and the 
results of the glueball spectrum. The study of the pseudoscalar channel using the topo- 
logical charge density operator will be discussed in Section 4, where we will also analyze 
the difference of the topological charge density operator from the conventional gluonic op- 
erators for the pseudoscalar glueball in previous quenched studies. Finally, we will give a 
summary and an outlook in Section 5. 


2 Lattice setup 
The gauge action we used is the tadpole improved gluonic action on anisotropic lattices [7]: 


Fo >> 5 Tr Pi; 1 Trkij 1 TrRji D9 4 ygT'r Poi 1 YgTr Rio 
ro a Oyu, 36 ygu 36 yguf ~ 19 u2 36 uf 


where Pj; is the usual plaquette variable and Rj; is the 2 x 1 Wilson loop on the lattice. 
The parameter us, which we take to be the forth root of the average spatial plaquette value, 
incorporates the usual tadpole improvement and y, designates the gauge aspect bare ratio 
of the anisotropic lattice, denoted as £o in our former quenched studies [20]. Although 
Yq suffers only small renormalization with the tadpole improvement [21], we have to tune 


it by determining the renormalized anisotropy ratio €,. As for the tadpole improvement 
parameter u; for temporal gauge links, we take the approximation up ~ 1 following the 
conventional treatment of the anisotropic lattice setup. 

We use the Wilson-loop ratios approach, with which the finite volume artifacts mostly 
cancel [22, 23]. We measure the ratios 


Wss (x, y) 


Rss = ~as Vs (yas) 2.1 

(x,y) maeta , (2.1) 
Wa(z, t) -as Vs (tat) 

Ra(z,t) = as Va (tar 2.2 

Ae Weta a 


and expect the spatial and temporal behaviors being the same at the correct €,. Therefore 
we find €, by minimizing 


L(q) = 5 (Reale, y) > Ralx, éau))" (2.3) 
ra (ARs) + (AR) 
where AR, and AR; are the statistical errors of Rss and Rst. We interpolate Rst(x, Egy) 
and its error with a cubic spine interpolation at non-integer égy. Since small x,y may 
introduce short-range lattice effects and large ones contribute only fluctuations, we scan 
and test different ranges and finally choose x,y € {2,3,4,5}. 
We adopt the anisotropic clover fermion action in the fermion sector [24]: 


p(z) (2.4) 


where Fv = tIm(P u (x)) and the dimensionless Wilson operator reads 


Wu = Va- sy 
Vule) = 5 [Uule) fle + w) — Uf — Fle p) 


Auf (x) = U l£) f(x + u) + Uh (@ — u) f(@ — u) — 2f (2). 


The bare fermion aspect ratio yf is also tuned to make sure that the measured aspect ratio 
Ef © Eg © € = 5. Ef is measured from the dispersion relation of the pseodoscalar and vector 
mesons 


pja? 
E(p)?a2 = ma? + 2 =. 
f 


where p = Ink, /Ls is the momentum on the lattice with periodic spatial boundary condi- 


(2.5) 


tions. 
The lattice spacings a, are set by calculating the static potential parameterized as 
V(r) = Vo +a/r+or. Using the Sommer scale pamameter ro’ = 410(20) MeV defined 


Table 2. Ground state meson spectrum 


Mpsat mps(MeV) my at my (MeV) Msat ms(MeV) 
0.07508(50)  650(4) | 0.114719) 993(16) | 0.1574(61) 1362(53) 
0.1119(4) 938(3) | 0.1388(12) 1164(10) | 0.1757(34) 1473(28) 


through r? we |r=r9 = 1.65, we can determine the ratio 


TO 1.65 +a 
= 2.6 
ds \ oa? (ae) 


where œ and ga? are derived from the fit to caluclated potential V(r) = V (fas) with f 


being the spatial distance in the lattice units. 

We generate two gauge ensembles on the 123 x 128 anisotropic lattice at 8 = 2.5 with 
the bare quark mass parameters mp = 0.05 and mp = 0.06. The pion masses on the two 
ensembles are measured to be 938 MeV and 650 MeV respectively. In the following, we 
use these m,,’s to label the gauge ensembles for convenience. The ensemble parameters are 
listed in Table 1. 

Apart from the pion masses, we also calculate the masses of the vector meson and 
scalar meson for calibration, which are listed in Table 2. We use the conventional J = 1 
vector and scalar quark bilinear operators as sink operators and the corresponding Gaussian 
smeared wall source operators to calculate the correlation functions. There is no ambiguity 
for the vector meson masses my’s since they are all below the two-pion threshold. For the 
scalar, we actually deal with ag whose two-body strong decay mode is mainly n'm (there 
is only one J = 0 pseudoscalar meson for Ny = 2, which is taken as the counterpart of the 
(approximately) flavor-singlet 1’ in the Ny = 3 case). At mz ~ 938 MeV, the calculated 
mass in ag channel is 1473(28) MeV, which must be the mass of ag since it lies below 
two-pion threshold and certainly below the 77 threshold. At Mme ~ 650 MeV, my is 
estimated to be my ~ 890 MeV (see below in Sec. 4), thus the mass value of 1362(53) 
MeV is also below the ņ'm threshold and can be taken as the mass of ag scalar at this 
pion mass. In order to calculate the J = 0 scalar meson mass, the disconnected diagrams 
(quark annihilation diagrams) should be considered. We have not done this yet, but as 
a raw assumption, we take the ag mass as an approximation to the mass of the isoscalar 
scalar meson. 


3 Numerical details 


In this work, the spectrum of the lowest-lying glueballs in three specific channels, namely 
scalar, tensor and pseudoscalar will be explored. The interpolating operators for these 
states are pure gluonic operators which have been extensively adopted in the previous 
quenched lattice studies. In other words, in each specific channel, no operators involving 


quark fields are included. This of course is only an approximation, assuming that the gluon- 
dominated state that we are after can be well-described by gluonic operators. Needless to 
say, mixing with the quark operators should be considered later on, especially for cases 
where the mixing is severe. For completeness, we briefly recapitulate the major ingredients 
of glueball spectrum computation in the following. One can resort to [9] for further details. 


3.1 Variational method 


The continuum SO(3) spatially rotational symmetry is broken into the discrete symmetry 
described by the octahedral point group O on the lattice, whose irreducible representations 
R are labelled as A,, A2, E, Tı, To, and have dimensionalities 1, 1, 2, 3, 3 respectively. 
Therefore, the lattice interpolation fields for a glueball of JPC quantum number should 
be denoted by RPC with R the irreducible representation of O which may include the 
components of J in the continuum limit. The parity P = + and the charge conjugation 


C = + can be realized by considering the transformation properties under the spatial 
reflection and time reversal operations. Since the octahedral group O is a subgroup of 
SU(2), the subduced representation of SU(2) with respect to O is reducible in general 
(for integer spin, this occurs for J > 2). Table 3 shows the reduction of the subduced 
representation of SU(2) up to J = 5. For instance, the scalar and pseudoscalar with J = 0 
states are represented by Aj, tensor states with J = 2 are reduced to direct sum of E and 
To, Le. (J = 2) 4 O = EQT. 


Table 3. Subduced representation of SU (2) with respect to O up to J = 5 


J 
R 0 12 345 
A, 1 0 0 0 1 0 
Ag 0 0 0 1 0 0 
E 0 0 1 O 1 1 
Tı 0 1 01 1 2 
To 0 01 1 1 1I 


As described in [9], we use Wilson loops (up to 8 gauge links) shown in Fig. 1. Each 
irrep R of group O can be realized by the specific linear combination of its 24 copies of a 
prototype Wilson loop under the 24 rotation operations of O. The combination coefficients 
of each R can be found in [9]. So each prototype may provide a different realization of R. 
On the other hand, the Wilson loops mentioned above can be built from smeared gauge 
links, such that different smearing schemes can provide more realizations of the gluonic 
operators. In practice, we have four different realization of each R by choosing different 
prototypes. For the smearing of gauge links, we adopt 6 smearing schemes by combining the 
single-link and double-link smearing procedures with different iteration sequences. Finally 
we have a set of 24 different gluonic operators, P, a = 1,2,...,24}, for each RPS. 


Paras ea ree 
AZ , 7 
— — m> m> 
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Figure 1. Prototypes of Wilson loops for the construction of glueball operators [8, 9]. 


Based on these operator sets, we use the variational method to get the optimized 
operators O(P) which mostly project to specific glueball states. In each symmetry channel 
R, we first calculate the 24 x 24 correlation matrix C® (t), 


CH et) = Sod Ere), a, 8 = 1,...,24 (3.1) 


T 


where pE is the vacuum-subtracted operator of of), 
Bt) = oH — ole (6)I0). (3.2) 


In practice, we only apply the vacuum subtraction to the operators in A[* channel. Sec- 
ondly, we solve the following generalized eigenvalue problem, 


CO (to)vi® = Ai(to)C™ Ove”, (3.3) 
where vif) is the i-th eigenvector, and A; = eo)! is the i-th eigenvalue where ™m,(to) 


is dependent on to and is close to the energy of the i-th state. For all the R channels, 
(R) 


Á gives the linearly combinational 
(R) 


i 


we use tọ = 1. It is expected that the eigenvector v 
coefficients of operators eo) to build an optimal operator P: ” which overlaps mostly to 


the i-th state, 


24 
E(t) = SoM GPW). (3.4) 
a=] 


3.2 Data analysis 
(R) 


In this work, the correlation function of the optimal operator ®; ” for the i-th state is 


calculated as 

(R R R 

Ct) = D oP t+ raa), (3.5) 
where we do the summation over the temporal direction to increase the statistics. Accord- 
ingly, the effective mass is defined as 


=(R) 
mË} (t) = In (so) (3.6) 
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(a) mz ~ 938 MeV (b) mz ~ 650 MeV 

Figure 2. Effective mass plateaus of a” (t) (red points) and ce (t) (blue points) in the R = Af* 
channel. The left and the right panel show the results at My ~ 938 MeV and m, ~ 650 MeV, 
respectively. The shaded bands are plotted with the best fit parameters using model of Eq. (3.9) 
in the illustrated time range. 


We divide the measurements into bins with each bin including 100 measurements. The 
statistical errors are obtained by the one-bin-eliminating jackknife analysis. 

For AL channel, the subtraction of the vacuum is very subtle. Even though we have 
O(10*) gauge configurations in each ensemble, when we perform the jackknife analysis 
above after subtracting the vacuum expectation values of the operator, we find there is still 
a residual (negative) constant term in the correlation function, which makes the effective 
mass M; e¢(t) going upward when t is large. This problem can be attributed to the large 
fluctuation of gauge configurations in the presence of sea quarks. To circumvent this 


difficulty, we adopt a vacuum-subtraction scheme by subtracting the correlation function 
C(t) with the shifted one C(t + ôt), 


— att ~ Att ~ Att 
C(t) = CAT) — EA" (tôt), (3.7) 
whose spectral expression is 
OAT a =D W; E emit) enmit =o it emit (3.8) 


ATT ; f ATT 
where W;;' is the spectral weight of the j-th state in C;* (t). 
constant term cancels with the spectrum unchanged. In practice, we take ðt = 5a,. 
We focus on the RPO = Aft, AJ", Ett, and T;* channels in this work. For all 


si a — Att 
these channels, the effective masses of CË(t) (red points) and CË(t) (blue points) ic i 
Aft 


Obviously, the possible 


for channel) are plotted in Fig. 2, 3, 4 and 5, respectively. In each figure, the left 
panel shows the result at m, ~ 938 MeV, and the right panel is for m, ~ 650 MeV. Even 
though we have a set of 24 operators for each channel, it is seen that the effective masses 
do not show plateaus from the very early time slices. This is very different from the case in 


the quenched approximation. One important reason for this is that, in each channel, the 


r 1st —— 
0.9 + 0.9 f+ 2nd 
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Figure 3. Effective mass plateaus of ae (t) (red points) and cy (t) (blue points) in the R = A7 t 
channel. The left and the right panel show the results at My ~ 938 MeV and m, ~ 650 MeV, 
respectively. The shaded bands are plotted with the best fit parameters using model of Eq. (3.9) 
in the illustrated time range. 
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(a) mr ~ 938 MeV (b) mz ~ 650 MeV 


Figure 4. Effective mass plateaus of Ge (t) (red points) and E®) (t) (blue points) in the R = E** 
channel. The left and the right panel show the results at My ~ 938 MeV and m, ~ 650 MeV, 
respectively. The shaded bands are plotted with the best fit parameters using model of Eq. (3.9) 
in the illustrated time range. 


spectrum of the full QCD is much more complicated than in the quenched approximation 
due to the sea quarks. This is true in principle, since qq states and multi-hadron states 
with the same quantum number do contribute to the corresponding correlation function in 
the presence of sea quarks. m 


Given the limited number of independent operators, our optimal operator ®; “ is 
actually not optimized as expected, namely, it does not only overlap to the i-th state but 


also to other states substantially. As seen in the effective mass plots, when mie), (t) tends 


to reach a plateau as t increases, m(t) decreases gradually and finally merges into this 
plateau at large t (within errors). Even though one can carry out the single exponential fit 
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Figure 5. Effective mass plateaus of C® (t) (red points) and oy (t) (blue points) in the R = T$" 
channel. The left and the right panel show the results at My ~ 938 MeV and m, ~ 650 MeV, 
respectively. The shaded bands are plotted with the best fit parameters using model of Eq. (3.9) 
in the illustrated time range. 


to the mass of the ground state in the plateau range roughly beyond t/a; ~ 6 or 7, the bad 
signal-to-noise ratio in this time range results in a large statistical error. Since we focus 
on the ground states in the present study, in order to get more precise result of the masses 
of the ground states, we adopt the following data-analysis strategy which also make use of 
the measured data in the short time range. In each channel, we carry out a correlated fit 
to a (t) and ana simultaneously through the following function forms, 


Gh (t) = Wi mit J WP emt, 
EL?) (t) = WEP emt WEP emt, (3.9) 


where the second mass term is introduced to take into account the contribution of the sec- 
ond state and higher states (of course, one can add more mass terms, but more parameters 
will ruin the data fitting due to the limited data points). In the fitting procedure, EW?) (t) 
and EL?) (t) are fitted in the same time window |tmin, tmax|. For each channel, we keep tmax 
fixed and varies tmin to check the stability and the quality of the fit. The fit results for the 
scalar (A{*), the pseudoscalar (A; *) and the tensor channels ( E++ and T" ) at the 
two pion masses are listed in Table 4 and 5. Except for tmin = 1 case in T t channel, 
all other fits are acceptable with reasonable x?/d.o.f. For all the four channels, the fitted 
parameters mı and Wj, are stable with respect to the varying tmin, while m2 decreases as 
tmin increasing gradually. This signals that our fitting model in Eq. (3.9) is not so good 
that we should include more mass terms to account for higher states, which, however, affect 
the second state more than the first state. Since we are interested only in the first states, 
we do not take mg seriously and treat it as an object accommodating the effect of higher 
states. We also illustrate the fit results in Fig. 2, 3, 4 and 5 with shaded bands. It shows 
that the fitting model Eq. (3.9) describes the data of the first states very well even at the 
very early time slices. For the second states, the fit model also fits the data more or less, 
especially in the small t region and the large t behavior. While in the middle t region, the 


Table 4. Fitted results using fit model Eq. (3.9) with different tmin at mz ~ 938 MeV. 


JPC Tiniti TM At m2 Qt Wi Wiz Wa) Wo2 x?/d.o.f 
At* 1 0.168(05) 0.545(07) 0.59(02) 0.37(02) 0.09(01) 0.86(01) 2.27 
2 0.166(07) 0.484(12) 0.58(03) 0.37(03) 0.05(01) 0.85(01) 0.50 
3 0.167(08) 0.471(19) 0.59(04) 0.35(05) 0.04(01) 0.83(02) 0.48 
4 0.167(11) 0.462(31) 0.59(06) 0.34(08) 0.04(02) 0.82(04) 0.52 
5  0.165(12) 0.508(52) 0.59(07) 0.42(12) 0.05(02) 0.94(14) 0.48 
Ay 1 0.306(08) 0.712(13) 0.65(02) 0.33(02) 0.13(2) 0.81(1) 2.24 
2 0.314(12) 0.650(24) 0.66(04) 0.28(05) 0.09(3) 0.80(2) 1.06 
3 0.302(18) 0.618(32) 0.60(07) 0.35(09) 0.07(03) 0.78(3) 1.09 
4 0.269(29) 0.530(47) 0.43(12) 0.54(13) 0.01(3) 0.69(4) 0.98 
Et++ 1 0,279(05) 0.695(09) 0.67(01) 0.31(01) 0.19(1) 0.77(1) 1.48 
2 0.280(06) 0.675(17) 0.66(02) 0.31(02) 0.18(2) 0.76(1) 1.48 
3 0.285(12) 0.582(32) 0.66(05) 0.27(06) 0.13(3) 0.71(2) 0.85 
4 0.281(18) 0.535(50) 0.63(10) 0.27(12) 0.09(5) 0.67(3) 0.86 
5  0.282(24) 0.508(75) 0.63(15) 0.25(20) 0.08(8) 0.64(6) 0.95 
Tt 1 0.282(04) 0.655(06) 0.63(01) 0.34(01) 0.151) 0.80(1) 4.47 
2 0.287(05) 0.597(09) 0.63(02) 0.32(02) 0.10(1) 0.80(1) 1.69 
3 0.300(08) 0.576(17) 0.70(04) 0.21(05) 0.10(2) 0.77(1) 1.08 
4 0.289(12) 0.563(31) 0.63(06) 0.31(08) 0.08(3) 0.77(3) 1.11 
5 0.275(17) 0.552(58) 0.55(09) 0.46(12) 0.07(4) 0.76(9) 1.13 


fitted results deviate somewhat from the data. This is understandable, since higher states, 
which do contribute, are missed in this model. This deviation actually contributes much 
to the x?. It is expected that the fitted mə is generally (much) higher than the mass of 
the second state. 

As shown in Table 4 and 5, most of the fits using different tmin are statistically 
acceptable and the masses of the first states are relatively stable. Therefore, for the final 
result of mı in each channel, we take tentatively the average value of mj ’s at different tmin 
weighted by their inversed squared errors. The statistical errors are accordingly derived. 
This averaging is illustrated in Fig. 6, where data points are the fitted result of mı at 
different tmin and the shaded bands are the averaged values with averaged errors. The 
results are also listed in Table 6. At the heavy pion mass m, ~ 938 MeV, m,(E*7) is very 
close to mı (1 t), as expected by the rotational symmetry restoration in the continuum 
limit where they correspond to the mass of the same 2** tensor state. However for the 


lighter m, ~ 650 MeV, the two masses deviate from each other by 200 MeV. Since the 


—10— 


Table 5. Fitted results using fit model Eq. (3.9) with different tmin at My ~ 650 MeV 


JPC Tiniti Mat m2 At Wii Wiz Wa) Wo2 x?/d.o.f 
Aj* 1 0.179(09) 0.544(12) 0.63(03) 0.35(03) 0.08(02) 0.86(02) 1.90 
2  0.159(14) 0.477(17) 0.54(05) 0.43(06) 0.04(2) 0.85(2) 0.94 
3 0.168(17) 0.49(30) 0.57(07) 0.39(08) 0.05(2) 0.86(03) 1.03 
4 0.178(20) 0.518(64) 0.63(09) 0.30(12) 0.06(4) 0.89(10) 1.16 
5  0.144(24) 0.527(91) 0.46(10) 0.75(20) 0.04(3) 1.00(27) 0.87 
A;* 1 0.291(08) 0.712(16) 0.61(02) 0.37(02) 0.12(2) 0.82(2) 2.04 
2  0.316(13) 0.689(33) 0.69(05) 0.25(06) 0.13(3) 0.79(2) 1.53 
3 0.333(19) 0.717(63) 0.77(09) 0.13(13) 0.16(5) 0.79(4) 1.65 
4  0.276(40) 0.580(86) 0.47(17) 0.56(18) 0.05(6) 0.70(7) 1.58 
E** 1 0.263(07) 0.629(08) 0.58(02) 0.39(02) 0.12(1) 0.82(1) 1.61 
2  0.266(10) 0.597(16) 0.58(03) 0.37(04)  0.10(2) 0.81(1) 1.30 
3. 0.255(13) 0.548(23) 0.52(05) 0.42(06)  0.06(2) 0.78(1) 1.08 
4 0.264(19) 0.546(40) 0.56(09) 0.36(11)  0.07(3) 0.77(5) 1.16 
5  0.296(45) 0.411(78) 0.77(48) -0.02(54) -0.12(28) 0.74(21) 0.74 
T5" 1 0.282(06) 0.671(08) 0.66(01) 0.31(02) 0.18(1) 0.78(1) 2.82 
2  0.301(10) 0.633(22) 0.71(04) 0.23(04) 0.16(3) 0.76(2) 1.07 
3 0.294(15) 0.588(32) 0.68(06) 0.25(08)  0.13(4) 0.75(2) 1.00 
4  0.287(24) 0.528(48) 0.63(12) 0.29(16)  0.07(6) 0.71(4) 1.05 
5  0.275(38) 0.513(76) 0.55(21) 0.40(28)  0.06(7) 0.71(8) 1.36 


Table 6. Final results for the masses of the lowest state we obtain in the Att, Ay +, Ett and 
T+ channels. These are the averaged values weighted by the inversed squared errors at each tmin- 


Ma(MeV) mi(Af*)(MeV) mi(E*+)(MeV) mi(T3°*)(MeV) mi(Ay*)(MeV) 
938 1397(25) 2342(25) 2392(25) 2559(50) 
650 1480(52) 2276(43) 2484(43) 2605(52) 


lattice spacings at the two pion massed are very close, the extent of the rotational symmetry 
breaking should be similar. We tentatively attribute this large deviation to the relatively 
small statistics at m, ~ 650 MeV, which is roughly one-half as large as that at m, ~ 938 
MeV (see Table 1). From Table 2 and Table 6 one can see that the masses of ground state 
scalar meson and our scalar glueball are very close to each other, this may indicate there 
are mixing between qq and the scalar glueball, which needs further investigation. 
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Figure 6. Fitted mı for At", AJ", Et+ and T;* are plotted with respect to tmin (the left panel 
for My ~ 938 MeV and the right panel for m, ~ 650 MeV). The values are expressed in the physical 
units inverted by the lattice spacing listed in Table 1. The color bands illustrate the averaged values 
weighted by the inversed squared errors at each tmin- 


3.3 Interpretation of the ground states 


Generally speaking, the two-point function of an interpolating operator O(t) with definite 
quantum numbers is usually parameterized as 


C(t) = (00H (0) |0) = $ 00n) (0t 0e, (3.10) 


n 


where {|n),n = 1,2,...} are eigenstates of Hamiltonian with eigen-energy Mn, which make 
up an orthogonal, normalized, and complete state set with 


So |n)(nl =1, Ha) = dn (3.11) 


For QCD on a Euclidean spacetime lattice, Mn take discretized values and the connection 
of these discretized energy levels to the relevant S-matrix parameters should be established 
through other theoretical formalisms, such as Ltischer’s. Here we would only focus on the 
physical meaning of the fitted masses of the lowest states. 

We take the scalar channel for instance. A hadron system of the bare states with 
the scalar quantum number JP? = 0*+ can be a bare scalar glueball |Go++), a bare qq 
scalar meson | fo), or even mr scattering states |nr). We simplify the matter further by 
assuming that the two adjacent states mix most, then we can only consider a two-state 
system composed of the ground state scalar glueball |G) and its adjacent state, which could 
be of nature |77) or | fo). This then yields the fitting model in Eq. (3.9) that we introduced 
previously. 

Despite the fact that glueball correlation functions in the unquenched QCD acquire 
more complicated spectrum decomposition than the quenched case, the mass of the bare 
glueball states |G’) can still be obtained by assuming the corresponding operators O couple 
weakly to other states. Therefore, it is naturally understood that the glueball spectrum 
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Table 7. We compare our results with previous results both from the quenched lattice QCD 
studies [8, 9] and the full-QCD study [13]. We average the masses of E**+ and T;‘* states to obtain 
the estimate of the 2** glueball mass. 


My (MeV) mo++ (MeV) mo++ (MeV) mo-+ (MeV) 


Ny =2 938 1397(25) 2367(35) 2559(50) 
650 1480(52) 2380(61) 2605(52) 
Ny =2+1 [13] 360 1795(60) 2620(50) = 
quenched [8] z 1710(50)(80) 2390(30)(120) 2560(35)(120) 
quenched [9] = 1730(50)(80) 2400(25)(120) 2590(40) (130) 


in our full-QCD lattice studies is similar to that in the quenched approximation. The 
difference is still visible, however, and it is most evident in the scalar channel where one 
would expect that this weak coupling assumption is not valid anymore. We compare the 
results in the present study with the previous quenched and unquenched results in Table 7. 


4 Further study on the pseudoscalar channel 


As presented in the last section, in the A;* channel, we obtain the bare pseudoscalar 
glueball mass m4-+ ~ 2.6 GeV at the two pion masses, which is compatible with the 
quenched results. A subtle question is that, due to the anomalous U,4(1) Ward identity, 


2 


On I$ (x) = 2mJ5(x) — T opie, (4.1) 
where g is the strong coupling constant, the gluonic operator €yypo F” FP? is usually called 
the topological charge density (up to a constant factor) and denoted by q(x), q(x) mixes 
with the flavor singlet pseudoscalar density J5(x), such that the gluonic pseudoscalar op- 
erator may have large overlap with the flavor singlet pseudoscalar meson (denoted by 1’). 
Since we are performing a Nf = 2 unquenched lattice calculation, one may expect that 
n’ can contribute to the pseudoscalar correlation functions. However, we do not observe 
the appearance of this state. So we would like to check the continuum form of the gluonic 
pseudoscalar (PS) operators involved in this work. Actually, in the construction of the 
gluonic pseudoscalar operators, only the spatially solid (instead of planar) Wilson loops 
(the last four prototypes in Fig. 1) are used, 


- -+ 
oa" (x,t) = So cat ReTr [Ro Wa(x, t) — PR o Wa(x, t)P7?], (4.2) 
REO 


where R stands for each rotation operation in O group, cAI is the combinational coefficients 
corresponding to the A, irreducible representation, Wa is any of the four prototypes made 
up of a specifically smeared gauge links. According to the non-Abelian Stokes theorem [25], 
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a rectangle Wilson loop Poe) of size a x b, with a,b small, can be expanded as 


1 
pe (x) = 1 + ab( F(x) + 5 (Du +bD,)Fu (£) 


i Z (20° D2 + 3abD Dy + 2b* D2) F u(x) 

4 5 (as + 2a°bD? D, + 2ab*D,,D? + b D3) Fuu (2)) 

+ (DPF) + 5 Fn) (Dy, + Dy) Fur (2) 

+ Fw) (a? D? + b?’ D?) Fiy(x)) 

+ = (ab Fi, + O((ab)*). (4.3) 


For simplicity, the factor ig is absorbed into the quantity F,,. The small ab expansion of 
Piy+v(x) is similar to Eq. 4.3 by replacing a and b with +a and b, respectively. Since 
the last four prototypes can be expressed as products of two rectangle Wilson loops, using 
the above relation one can obtain the leading term of the pseudoscalar operator, 


Alt 
Qa’ (x,t) x eijk TrBi(x, t)D;Bk(x, t) + O(a?), (4.4) 
which is obviously different from the anomaly part of the anomalous WI, epvpo FH” (x) FP? (x) œ 
E(z)- B(x). This may imply that the two operators couple differently to specific states. In 
order to clarify this, we would like to use q(x) as the operator for pseudoscalars and study 
their spectrum. The correlation function of q(x) is 


2.0E-4 = 2.0E-4 =r 
1.0E-4 ae 1.0E-4 Eos 
0.0E0 ines 0.0E0 peer 
S  -1.0E-4 5 ~1.0E-4 
-2.0E-4 -2.0E-4 
-3.0E-4 -3.0E-4 
-4.0E-4 -4.0E-4 
(a) Mr ~ 938 MeV (b) mz ~ 650 MeV 


Figure 7. The correlation function C,(r) of topological charge density in terms of the four dimen- 
sional Euclidean distance (the left panel for Mmr ~ 938 MeV and the right panel for my ~ 650 MeV. 
Different curves correspond to C,(r) at different Wilson flow time t = 0.2, 0.3, 0.4 and 0.8. 


Calz — y) = (a()a(y))- (4.5) 
Since the topological susceptibility X+ = vA f d*xd*yC,(x—y) > 0 (V4 is the four-dimensional 
volume of the Euclidean spacetime), and Cy(a — y) < 0 for r = ||x — y|| > 0, intuitively 
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C(x — y) can be expressed as 
Cq(a — y) = Ad*(a — y) + Cy(a — 9), (4.6) 


where C,(x — y) is negative for r > 0. On the Euclidean spacetime lattice with a finite 
lattice spacing, the delta function will show up a positive kernel with a width of a few lattice 
spacings, and C,(x — y) has a negative tail contributed from C,(x — y). Because q(x) is 
a flavor singlet pseudoscalar operator, it is expected that C,(2 — y) would be dominated 
by the contribution of the lowest pseudoscalar meson in the large r range and can be 
parameterized as [26] 


Cy(r) = NTE Ki (mpsr), (4.7) 


where N is an irrelevant normalization factor, mpg is the mass of the lowest pseudoscalar, 
and K(z) is the modified Bessel function of second kind, whose asymptotic form at large 


T 3 3 
Kiı(z) ~ ia “(1+ an” larg z| < 57 (4.8) 


Therefore, one can obtain mpg by fit the negative tail of C(x — y) in the large r range using 


|z| is 


the above functional form. This has been actually done by several lattice studies in both 
the quenched approximation [19] and full QCD calculations [18]. In the quenched approx- 
imation, the extracted mpg = 2563(34) MeV is in good agreement with the pseudoscalar 
glueball mass mps = 2560(35) MeV. This is what should be, since the hadronic excitations 
of a pure gauge theory are only glueballs. In the full-QCD study with Nf = 2+ 1 and pion 
masses close to the physical my, mpg is obtained to be 1013(117) MeV, which is consistent 
with the mass of the physical 1’. 

In this work, we adopt the similar strategy to that in [18]. The topological charge 
density q(x) is defined by the spatial and temporal Wilson loops (plaquettes) as conven- 
tionally done. We use the gradient Wilson flow method as smearing schemes to optimize 
the behaviour of topological charge density correlator [18, 27]. The Wilson flow provides a 
reference energy scale Ja [28]. In practice, we use the code published by the BMW collab- 
oration [29] to evaluate the topological charge density. Fig. 7 shows the topological charge 
density correlator for My ~ 938 MeV and m, ~ 650 MeV at flow time t = 0.2, 0.3, 0.4, 0.8 
respectively. On our lattices, these t values correspond to v8t ~ 0.15,0.18,0.21 and 
0.30 fm. As shown in the figures, at large flow time, the correlator is almost positive 
which implies that the gauge fields are over smeared. We then use the correlation function 
at flow time t = 0.4, which seems to have better signal-to-noise ratio and is smoother in 
our fitting range, to extract the mps. 

We fit the negative tail region of C (r) with Eq. 4.7, and the fitting range is determined 
by defining the effective mass mer ¢(r) of n through 


C(r+Ar)_ or Kilmess(r + Ar)) 
C(r) r+Ar Ky (meg fr) 


(4.9) 


where we take Ar ~ 0.05 fm and C(r) is averaged in [r,r + Ar], then mer f(r) is obtained 
by solving Eq. 4.9 numerically. Fig. 8 shows the effective mass plateaus for flow time 
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Figure 8. The effective mass plateaus for Wilson flow time t = 0.2,0.3,0.4 at Mr ~ 938 MeV (left 
panel) and my ~ 650 MeV (right panel), respectively. The color bands show the fitted masses with 
errors and also the fit range. 


t = 0.2,0.3, 0.4 whose errors are estimated by jackknife method. Due to the small lattices 
and the periodic boundary condition we are using, we suffer from serious finite volume 
effects in fitting mpg since the positive kernel extends about 2a, and the fitting range 
should be far enough from L/2 = 6asg. 

Table 8 shows the fitted result of mpg at the two pion masses, which show explicitly 
the quark mass dependence of mpg. For the large uncertainties, we do not take the fitted 
values of mpg too seriously but obtain the impression that they are far below the pseu- 
doscalar glueball mass in the quenched approximation and that in this study. Anyway, the 
mpg’s seem reasonable to some extent. Physically, the large 7 mass is acquired through 
the interaction of sea quark loops and the mechanism can be described by the following 
expression, 


1 1 1 1 1 
= (aie +m? m? te (4.10) 
q? — m2, a ee a i 
such that 
ma sm? + må, (4.11) 


where me is related to the topological susceptibility x+ through Witten- Veneziano mecha- 
nism [30, 31] as 


4N 
mô = a Xo (4.12) 
where fr is the decay constant of m. For our case of Ny = 2, if we take the values 


Xt = (180 MeV)*, fr ~ 150 MeV for m, ~ 650 MeV and fp ~ 200 MeV for m, ~ 938 
MeV, m2 is estimated to be approximately (610 MeV)? and (460 MeV)?, respectively. Thus 
the y’ mass can be derived as my ~ 890 MeV for mz ~ 650 MeV, and my ~ 1045 MeV 
for Mmr ~ 938 MeV. These values are not far from the mps’s we obtained. 

Even though these are very coarse calculations, from which we can see that there does 
exist in the spectrum a flavor singlet gq pseudoscalar meson corresponding to the 7 meson 
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Table 8. The fitting details for ņ meson mass from topological charge density correlator at 
Mz ~ 938 MeV and my, ~ 650 MeV respectively. 


Mr fit range(as) Myas | My (MeV) | x?/dof 
938 MeV | 3.74-5.92 | 0.856(21) | 1481(36) | 1.01 
650 MeV | 3.87-5.48 | 0.514(22) | 890(38) | 1.43 


in the real world, which can be accessed by the topological charge density operator. Along 
with the observation in the calculation of glueball spectrum, this proves to some extent 
that our operator for the pseudoscalar glueball couples very weakly to the qq meson state 
but almost exclusively to the glueball states. 


5 Summary and conclusions 


The spectrum of the lowest-lying glueballs are investigated on the lattice with two flavors 
of degenerate quarks. We generate large ensembles of gauge configurations on anisotropic 
lattices at two pion masses, Mr ~ 650 MeV and Mmr ~ 938 MeV. We focus on the ground 
states of the scalar, pseudoscalar and tensor glueballs, which are measured by gluonic 
operators constructed from different prototypes of Wilson loops. The variational method 
are applied to obtain the optimal operators which couple dominantly to the ground state 
glueballs. 

In the tensor channel, we obtain the ground state mass to be 2.367(35) GeV and 
2.380(61) GeV at m, ~ 938 MeV and 650 MeV, respectively. In the pseudoscalar channel, 
when we use the gluonic operator whose continuum limit has the form of ¢€;;,.7'rB;Dj Bz, 
we obtain the ground state mass to be 2.559(50) GeV and 2.605(52) GeV at the two pion 
masses. These results are in agreement with the corresponding glueball masses in the 
quenched approximation and show little dependence of pion masses. In the scalar channel, 
the ground state masses extracted from the correlation functions of gluonic operators are 
determined to be around 1.4-1.5 GeV, which are close to the ground state masses from the 
correlation functions of the quark bilinear operators. In all the channels considered, the 
spectrum is similar to the quenched case and the mixing between glueballs and conventional 
mesons will be explored in the future. 

We also investigate the pseudoscalar channel using the topological charge density as 
the interpolation field operator, which are defined through Wilson loops and smeared by 
the Wilson flow technique. The masses of the lowest state derived in this way are much 
lighter (around 1GeV) and compatible with the expected masses of the flavor singlet qq 
meson. This provides a strong hint that the operator ¢€;;,7'rB;Dj;B, and the topological 
charge density (proportional to TrE - B) couple very differently to the glueball states and 
qq mesons. 

Admittedly the lattice volume we used is relatively small and the continuum limit re- 
mains to be taken, our current results are still helpful to clarify some aspects of unquenched 
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effects of glueballs and serves as a starting point for further future studies. 
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